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1 .  Summary 


The  present  work  developed  a  three-dimensional  composite  velocity 
solution  procedure  for  the  reduced  Navier- Stokes  equations.  The  scheme 
developed  is  consistent  with  both  asymptotic  boundary  layer  and  inviscid 
flow  theories.  The  composite  velocity  scheme  has  been  differenced  to 
provide  second  order  accurate  solutions.  Global  storage  is  required  for 
only  three  variables.  If  storage  becomes  a  limiting  factor,  this  may  be 
reduced  to  one  variable  by  coupling  the  r)-momentum  and  energy  equations 
to  the  rest  of  the  system.  In  separated  flow  regions,  we  found  that 
only  the  derivatives  needed  to  be  upwinded. 

To  solve  the  resulting  finite  difference  equations,  a  standard  plane 
relaxation  procedure  has  been  used  with  a  consistent  coupled  strongly 
implicit  procedure  to  invert  the  system  at  each  cross  plane.  The  CCSIP 
developed  is  an  approximate  LU  decomposition  that  provides  a  solution  at 
each  cross  plane  that  is  second  order  accurate  with  no  iterations  on  the 
inversion  scheme  required.  This  scheme,  coupled  with  a  second  order 
accurate  quasi-linearization  of  nonlinear  terms  provides  a  global 
solution  procedure  that  is  second  order  accurate  and  requires  no  local 
iterations  at  each  cross  plane. 

Results  for  transonic/turbulent  flow  over  both  axisymmetric  and 
nonaxi symmetric  geometries  have  been  presented.  The  composite  velocity 
scheme  is  an  efficient  method  for  solving  many  three-dimensional 
problems.  Several  possible  improvements  to  the  algorithm  have  been 
pointed  out. 
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2.  Introduction 


As  the  field  of  computational  fluid  dynamics  has  matured,  the  need 
for  solution  procedures  that  can  efficiently  solve  three-dimensional 
flows  with  complex  viscous/inviscid  interactions  has  become  apparent. 
One  technique  for  solving  these  problems  is  to  obtain  solutions  to  the 
time -dependent  Navier- dtokes  equations  using  the  time  asymptotic  limit 
for  the  steady-state  solution.  While  the  rapid  advances  in  computer 
technology  have  made  this  type  of  solution  procedure  possible,  they  are 
still  costly  in  both  computer  time  and  resources. 

The  reduced  Navier-Stok.es  (RNS)  equations  provide  an  attractive 
alternative  to  the  full  time -dependent  Navier-Stokes  equations  for 
viscous/inviscid  interacting  flows.  This  set  of  equations  has  been 
developed  to  simulate  large  Reynolds  number  (Re)  asymptotic  behavior 
while  retaining  a  single  composite  set  of  equations.  The  RNS  equations 
are  a  composite  of  the  full  Euler  and  second-order  boundary  layer 
systems.  The  terms  neglected  in  the  Navier-Stokes  equations  are  higher 
order  in  Re  in  an  appropriate  'streamline'  coordinate  system. 

The  asymptotic  nature  of  the  RNS  equations  may  be  further  exploited 
by  the  use  of  a  composite  representation  for  the  flow  velocities.  In 
the  composite  velocity  formulation,  a  combination  of  a  multiplicative 
and  additive  composite  is  defined  In  the  spirit  of  matched  asymptotic 
expansions.  Viscous  or  rotational  velocities  U  and  W  and  a  pseudo¬ 
potential  or  inviscid  velocity  (<J>x>  4>^,  4>^)  are  specified.  The 

composite  velocity  formulation  is  thus  consistent  with  both  asymptotic 
inviscid 
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flow  and  boundary  layer  theory.  For  inviscid  flows,  the  continuity 
equation  reduces  to  the  full  potential  equation,  and  the  momentum 
equations  are  identically  satisfied.  In  viscous  flow  regions,  the 
velocities  are  determined  from  the  continuity  equation  and  the  axial  and 
crossflow  momentum  equations.  The  normal  momentum  equation  provides  a 
viscous  total  pressure  correction.  The  composite  velocity  formulation 
has  been  applied  previously  for  the  two-dimensional  Navier-Stokes 
equations  for  both  incompressible  [1]  and  compressible  flows  [2-4].  In 
previous  work  by  the  present  author,  the  composite  velocity  formulation 
was  applied  to  the  two-dimensional  Euler  and  RNS  equations  [5-6]. 

The  splitting  of  the  velocities  into  a  composite  of  viscous  and 
inviscid  velocities  provides  several  advantages  over  primitive  variable 
formulations.  Since  the  full  potential  equation  is  the  basic  equation 
for  the  inviscid  flow  in  the  composite  procedure,  the  extensive 
knowledge  base  for  solving  this  equation  for  transonic  flows  may  be 
exploited  when  developing  the  solution  procedure  for  the  composite 
velocity  system.  In  particular,  techniques  for  transonic  flow  developed 
for  the  full  potential  equation  may  be  applied  to  the  composite  velocity 
form  of  the  RNS  equations.  In  the  present  work,  the  Enquist-Osher  flux 
biasing  procedure  [7]  is  used  in  supersonic  regions.  The  ability  to 
identify  viscous  and  inviscid  velocities  allows  these  terms  to  be 
treated  numerically  appropriate  to  the  asymptotic  character  of  the 
terms,  i.e.,  marching  of  the  viscous  velocities  representing  the 
parabolic  character  of  the  boundary  layer  equations  and  central 
differencing  for  the  inviscid  velocities  representing  the  elliptic 
nature  of  inviscid  subsonic  flow.  For  separated  flows,  the  two- 
dimensional  composite  velocity  formulation  does  not  require  upwinding  in 
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separated  region,  and  the  three "dimensional  formulation  requires 


.ie 

upwinding  for  only  the  azimuthal  W 


f 


convective  derivative. 


The  composite  velocity  procedure  is  applied  to  the  solution  of 
nonaxi symmetric  afterbody  flows.  Solutions  to  these  types  of  problems 
can  be  very  important  in  aft-fuselage -engine  integration  for  fighter 
aircraft  design.  The  complex  viscous/inviscid  interactions  in  these 
flows  provide  a  computational  challenge  to  any  flow  solver.  The 
composite  velocity  procedure  provides  an  efficient  alternative  to  some 
of  the  Navier- Stokes  procedures  currently  being  used  to  solve  these 
problems  [8-9]. 


4 


3.  The  Governing  Equations 


The  governing  equations  to  be  solved  are  given  by  the  three- 
dimensional,  steady,  compressible,  RNS  approximation.  These  equations 
have  been  formulated  for  a  general,  nonorthogonal ,  curvilinear  system  of 
coordinates.  This  allows  the  use  of  appropriate  'streamline'  coordinate 
systems  and  places  no  restrictions  on  the  type  of  grid  to  be  used. 

The  governing  equations  are  now  given  in  their  composite  velocity 
formulation.  In  the  spirit  of  matched  asymptotic  expansions,  the  flow 
velocities  are  written  as 


u  =  (U+l)(g‘  +  g'1*  +g‘ 1<J>^)=(U+l)ue 

v  =  (g2'^  +  g"*  +  g”^) 


w  =  (U+l)(g3  l<t>£+g32<i>T+g33<i^)  +  W 


=  (U+l)we+  U 


The  composite  representations  of  u  and  w,  the  axial  and 
velocity  components,  contain  two  types  of  terms,  irrotational 
potential  velocities,  e.g.  ,  <J>^,  and  viscous  velocities  U 


(la) 

(lb) 

(lc) 


crossflow 
'pseudo' 
and  W. 


Since  the  change  in  v  across  the  boundary  layer  is  small,  the  normal 
velocity  is  determined  solely  by  the  'pseudo'  potential.  This 
particular  composite  form  has  the  advantage  that  no  additional  unknowns 
are  introduced,  and  the  boundary  conditions  remain  easy  to  implement. 

By  substituting  equations  l(a-c)  into  the  RNS  equations,  the 
following  system  is  obtained: 
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Continuity 

(pv'g(U+l)ue)  +  |^(pv/gv)  +  g(U+l)we)  +  ^(pv^gW)  =  0  (2) 


f -Momentum  Equation 


lfj(P^g(U’+U)u’)  +  ^(p^gUvu,)*  f7(p^g(UJ+U)uewe)  (3) 


P/ g 


V  3rjv^6  e 


e  e 


+  ^-(p^gUWu  )]  +  Uu  u  +  Uw  u  +  Wu  +  curvature  terms  = 
at  e  e  ef  6  e£ 


-  g"[Gf-TS,]  -  g'*[G  -TSJ  -  g '  *  [G  -TS  ]  +  -*z  (viscous  terms) 

f  f  V  V  <  <  pV g 

f -Momentum  Equation 


~  [|r(Wg(U2+U)uewe)  +  gT(P'/g(U+l)Wue)  +  f-(p^gUvwe) 
p'/g  4 

+  f^(p^gvW)  +  |^(pVl(U2+U)w*)  +  “(/*/£(  2U+l)Wve) 


(4) 


+  f~(pVgW  )]  +  Uu  w  +  Uw  v  +  Ww  +  curvature  terms  = 

PC  e  e^  e  ef  ef 

-g3‘[G^-TS?]  -  g**[GTj-TS7j]  -  g,3[C^-TS^]  +  -±z  (viscous  terms) 


g 


Normal  Momentum  Equation 


Uu  v.  +  (Uw  +W)vt  +  curvature  terms  = 

e  (  e  C 


■g*  * [G^-TS^] 


(5) 


z”lCVmTSV]  ’  s’^V^C3 


where 


v  p  .  1  ij  i  j 
=  ~Zr  *  +  ~  g  JU  UJ  . 
y-1  p  2  °  e  e 
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The  energy  equation  may  also  be  included  in  the  composite  velocity 
system  as  an  equation  for  temperature.  In  the  present  calculations, 
however,  it  is  replaced  with  the  assumption  of  constant  total  enthalpy. 
This  is  a  valid  approximation  for  the  Mach  number  regime  being 
considered.  The  new  variable,  G,  that  appears  in  these  equations  is 
similar  to  the  total  (or  Bernoulli-like)  pressure.  G  is  not,  however, 
assumed  to  be  constant  but  is  calculated  from  the  solution  procedure. 

In  viscous  regions,  the  composite  velocity  formulation  is 
representative  of  the  reduced  form  of  the  NS  equations.  The  continuity, 
£ -momentum,  and  (--momentum  equations  determine  U,  W,  and  <I>.  The 
temperature,  T,  is  obtained  from  the  energy  equation  (or  total  enthalpy 
H)  and  the  total  pressure  correction  is  determined  from  the  rj-momentum 
equation.  In  the  inviscid  limit  as  U,W-0,  the  continuity  equation 
reduces  to  the  full  potential  equation  and  the  Bernoulli  relation, 
G=constant,  is  recovered  from  the  momentum  equations.  Thus  equations 
(3-5)  are  identically  satisfied,  and  the  composite  velocity  system 
reduces  to  the  expected  representation  for  an  inviscid,  irrotational 
flow. 

The  composite  velocity  scheme  is  formulated  to  provide  full  Euler 
solutions  for  the  outer  inviscid  flow.  If  the  momentum  equations  are 
solved  in  the  present  nonconservation  form  (3-4),  however,  the  full 
potential  solution  is  recovered  for  the  outer  inviscid  flow  instead; 
i.e.,  neither  entropy  nor  vorticity  is  generated  in  the  outer  inviscid 
flow.  This  includes  flows  with  captured  shock  waves.  To  capture  the 
rotational  or  Euler  shock  wave  in  transonic  flows,  the  (-momentum  and  (- 
momentum  equations  must  be  written  in  a  quasi -conservation  form  [6]  as 
follows: 
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A -Momentum  Equation 

~~  [|t(/VI(U2+2U)u^)  +  |-(p^iUvue)  +  f-(/Vl(U2+2U)uewe)  (6) 

g 

+  “(pyg(U+l)Wue)]  +  curvature  terms  =  -  ^[g''p^  +  g'^  +  g1^] 

-  “ ~  [|j(P^l gu2)  +  ^(pVgvu,)  +  |^(pv'guewe)]  -  curvature  terms 

+  — ^  {viscous  terms) 
pJ  g 


/•-Momentum  Equation 

“  [|T(p^i(U2+2U)ueve)  +  |r(pv'g(U+l)Wue)  +  §-(p^Uvwe)  +§-(p</gvW) 

ps/g  ?  * 

+f^:(P^g(U2+2U)v2)  +  |^;(2p^I(U+l)Wve)  +  ^(p^W2)]  (7) 

1,11  12  IS, 

+  curvature  terms  =  -  ~[g  +  g  p^  +  g  p^J 

2  ^  __  ^  ^  j 

-  —I  [^j(p^guewe)  +  g^(py gvwe)  +  ^(Wgwe)]  *  curvature  terms 

+  — ^  {viscous  terms) 

pJ g 

The  rj-momentum  equation  (5)  is  retained  in  nonconservation  form  and 
represents  the  Crocco  equation  that  relates  entropy  to  vorticity. 

With  the  equations  in  'conservation'  form,  the  correct  entropy  rise 
at  the  shock  wave  will  now  be  generated;  however,  spurious  entropy  is 
also  generated  in  nonshock  regions.  In  the  present  technique,  a  simple 
solution  to  this  problem  is  available.  The  nonconservative  form  of  the 
momentum  equations  (3-5)  produces  no  entropy,  but  accurately  convects 


the  entropy  or  vorticity  present  in  the  flow.  Therefore,  this  form  of 
the  axial  momentum  equation  is  used  everywhere  except  in  the  shock 
region  where  equations  (6,7)  are  required.  This  leads  to  a  solution 
procedure  with  the  desirable  feature  of  generating  the  correct  entropy 
rise  at  the  shock  wave  but  not  creating  spurious  entropy  in  other 
regions  of  the  flow.  The  advantage  of  nonconservative  equations  away 
from  shock  waves  and  combined  conservative/nonconservative  systems  has 
been  discussed  in  several  studies  by  other  investigators  [11]. 
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4.  Boundary  Conditions 


Boundary  conditions  for  the  composite  velocity  formulation  may  be 
easily  specified.  At  solid  surfaces,  the  no-slip  condition  and  zero 
injection  conditions  are  used  giving  U=-l,  W=0,  and  v=0.  The  upper 
boundary  is  to  be  sufficiently  far  from  the  body  that  the  flowfield  is 
undisturbed.  This  yields  freestream  boundary  conditions  U=0,  W=0, 
<i>=<t)£s,  T=l,  S=0.  At  the  outflow,  ,  only  one  boundary  condition 

needs  to  be  applied,  4>^=0.  This,  in  effect,  assumes  a  weak 

viscous/inviscid  interaction.  At  the  inflow  boundary,  uniform  flow  is 
assumed  for  all  values  except  at  solid  boundaries  where  the  no-slip 
condition  is  applied.  For  all  geometries  presented  in  this  paper,  two 
planes  of  symmetry  exist  in  the  azimuthal  direction.  In  each  of  these 
planes,  symmetry  conditions  are  applied  for  U,  4>,  T,  S,  with 
antisymmetry  conditions  for  W. 
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5.  Numerical  Procedure 


The  numerical  procedure  developed  for  the  composite  velocity 
equations  is  motivated  by  numerical  procedures  for  potential  flows,  <t>, 
and  for  boundary  layer  flows,  U,  W.  In  the  composite  equations,  the 
upstream  or  acoustic  influence  is  contained  wholly  in  the  potential  part 


of  the  equations,  i.e.,  4> 


A  standard  plane  relaxation  procedure 


similar  to  that  for  solutions  to  the  full  potential  equation  is  used  for 
all  $  derivatives.  The  viscous  velocities,  U  and  W,  are  treated  as 
boundary  layer  velocities  and  may  be  marched  except  in  regions  of 
separated  flow.  The  continuity,  f -momentum,  and  { -momentum  equations 
are  solved  as  a  coupled  set  of  equations  for  U,  W,  and  A<1>.  The  t)- 
momentum  equation  and  energy  equation  are  solved  uncoupled  from  the 
system  for  entropy,  S,  and  temperature,  T. 


5.1  Finite  Difference  Equations 


The  composite  velocity  equations  are  differenced  so  that  second 
order  accuracy  is  obtained  for  all  terms.  First  order  accurate 
differencing  for  the  streamwise  derivatives,  and  W^.,  is  left  as  an 

option.  The  difference  form  of  the  continuity,  ^-momentum,  and  (- 
momentum  equations  is  obtained  at  (i,j,k),  with  central  differences  at 
half  points.  The  derivatives  of  4>  are  then  central  differenced  at  the 
half  point  locations.  This  provides  the  usual  three  point  central 
difference  for  <J> _ ,  and  <t>  .  In  the  n-/-  cross  plane,  the  values  of 

U  and  W  at  the  half-points  are  determined  by  averaging  the  values  at 
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neighboring  grid  points.  In  order  to  provide  the  proper  upwind  bias, 
consistent  with  the  boundary  layer  marching  character  of  U  and  W,  the 
values  of  U  and  W  at  i+1/2  and  i-1/2  are  represented  with  upwind 
approximations  as  follows: 

Ui+l/2  '  “l  +  2  <8a> 

“i-1/2  ‘  “i-i  +  2  e“i-l<Ui-l-  “i-2>  <8b) 

where  e=0  provides  a  first  order  accurate  representation  and  e=l 
provides  a  second  order  accurate  representation. 

In  separated  flows,  the  representation  of  anc*  ^i-1/2  must 


modified  to  provide  the  proper  upwinding  for  the  derivative.  This  is 


required  since  W  appears  as  an  additive  term  in  the  composite 
representation  for  w,  rather  than  a  term  multiplying  an  inviscid 
velocity,  e.g. ,  U.  Therefore  the  only  upstream  influence  in  separated 
regions  comes  through  upwinding  of  the  derivative.  The 


representations  for  anc^  ^i-1/2  become 

“i+1/2  '  SM*'"i  +  2eW  "i-i”  +  SP*I"i+l  +  2^  <Ui+l-  Ui+2>J  <9a> 

"i-1/2  ’  SM*t"i-l  +  K-l("i-l-  "i-2,l  +  SP*t"i  +  2^  <"i-“i+l>l  <9b> 


where 


SM  =  |[2  +  SGN(U1 _x+l)  +  SGN(U^+1)] 
SP  =  |[2  -  SGNCU^j+l)  -  SGN(U1+1)]. 


12 


In  equations  (9a, b)  a  backward  difference  for  results  if  SM=1  and 


SP=0,  and  a  forward  difference  results  when  SM=0  and  SP=1.  At 
separation  and  reattachment  points,  SM=l/2  and  SP=l/2,  and  a  central 
difference  is  obtained  for  W^.  This  provides  for  a  smooth  transition 


between  backward  and  forward  differencing  at  separation  and  reattachment 
points. 

The  first  order  boundary  layer  viscous  terms  are  treated  implicitly 
and  have  been  linearized  by  assuming  all  inviscid  velocities,  ue>  v,  and 


wg  known.  Since  the  second  order  boundary  layer  terms  are  small ,  these 
terms  have  been  included  as  an  explicit  correction,  with  the  inviscid 


velocities  calculated  from  the  previous  global  iteration  and  U  and  W 


determined  from  previous  marching  steps.  The  combinations  G^-TS^.,  G^- 


and  G^.-TS^. 


also  have  a  minimal  influence  on  the  solution  of  the 


governing  equations.  These  terms  have  been  uncoupled  from  the  momentum 
equations  and  are  assumed  known  from  the  previous  global  iteration. 

Certain  terms  must  also  be  lagged  so  that  a  stable  relaxation 
procedure  is  obtained.  In  the  global  relaxation  process,  the  only  term 
that  introduces  <l>  at  two  different  iteration  levels  is  <J>^.  All  first 


derivatives  of  4>  In  the  streamwise  direction, 


and  the  cross 


derivatives  3>. 

to 


are  taken  from  the  previous  global  iteration. 


All 


r)-{  cross  derivatives  are  assumed  known  so  that  a  five  point 
computational  star  is  maintained  in  the  cross  plane. 

When  solving  the  difference  form  of  the  equations,  errors  arise  in 
the  freestream  because  of  incomplete  cancellation  of  metric  terms[12]. 
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Essentially,  a  group  of  terms  that  analytically  should  be  zero 
introduces  an  error  owing  to  the  differencing  of  metric  terms.  The 
following  residual  term  represents  the  freestream  error  for  the 
continuity  (i.e.,  full  potential)  equation 


FSR  =  fj(^efs)  +f^&vfs>  +  <10> 

where  the  freestream  velocities  are  calculated  from  ^fs*  To  exactly 

reproduce  freestream  results,  this  error  term  is  subtracted  cut  of  each 
equation. 

The  ^-momentum  equation  is  solved  uncoupled.  It  is  treated  as  a 
first  order  differential  equation  for  entropy,  S.  The  difference  form 
of  the  equation  is  centered  at  (i-1/2,  j+1/2)  with  all  terms  being 
central  differenced  at  this  location.  The  terms  arising  from 


nonorthogonal  grids,  G^-TS^  and  G^-TS^,  are  taken  from  the  previous 

global  iteration.  All  other  terms  in  the  equation  involve  U,  W,  4>, 
which  have  already  been  calculated. 

Employing  the  assumptions  for  the  plane  relaxation  procedure,  the 
final  3x3  system  of  difference  equations  to  be  solved  at  each  cross 
plane  are  written  as 

Af^+1.  ,  .  +  bV?+1.  +  +  D'P?+1.  ,  .+  ET?+1,  . 

i,j-l,k  i,j,k  i , j+1 , k  i,j,k-l  i,j,k+l 


G  -  -  Hr"  .  , 

i-l,j,k  i+1, j,k 


(ID 


where  n+1  represents  the  current  iteration  level;  n  represents  the 


previous  global  iteration,  and  ¥=[U, W,  A<i>]  .  Once  equation  (11)  has  been 


solved,  the  values  of  j  k  ^i-1/2  j  k  are  t^en  updated  from 

the  energy  and  T)-momentum  equations,  respectively. 


5.2  Consistent  Coupled  Strongly  Implicit  Procedure 


To  solve  the  system  of  difference  equations  (11)  in  an  efficient  and 
robust  manner,  a  solution  algorithm  that  is  noniterative, 
unconditionally  stable,  and  spatially  consistent  is  required.  A 
noniterative  algorithm  is  important  since  this  algorithm  is  applied  at 
each  cross  plane  in  the  marching  direction  and  a  number  of  global 
iterations  is  required  to  obtain  a  converged  solution.  An 
unconditionally  stable  and  spatially  consistent  scheme  is  desirable 
since  grid  stretching  in  the  marching  direction  allows  for  greater 
flexibility  in  resolving  important  flow  features  with  a  minimum  number 
of  grid  points.  A  limitation  on  the  step  size  in  the  marching  direction 
owing  to  stability  restrictions  or  the  need  to  maintain  the  spatial 
consistency  of  the  scheme  would  make  the  number  of  planes  required  in 
this  direction  prohibitive. 

A  solution  algorithm  that  satisfies  the  above  requirements  is  the 
consistent  coupled  strongly  implicit  procedure  (CCSIP)  developed  by 
Khosla  and  Rubin[13-14] .  The  consistent  CSIP  algorithm  may  be  derived 
from  the  following  recursion  relationship: 


^i.j.k  ^j.k 


+  ET.  .4'.  .  1  .  +  FT .  ,  ¥ .  .  , 

j,k  i,j-l,k  J,k  j,k+l 


where  GM.  ,  is  a  column  matrix  and  ET.  ,  and  FT.  , 

J » K  J  >  *  J  >  k 


matrices. 


Equation  (12)  is  substituted  into  equation  (11)  for  4',  .  , 

1  t  J  +  K 
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'J'i  .  ^ The  corner  points  introduced  in  this  substitution  are 


represented  as  follows 
,.n+l  * 


.  n*r  1  _  *  x  ^  x 

fi,j+l,k+l  =  Vi,j+l,k+l"  “^i.j+l.k  +  *i, j.k+l'  ^i, j,k} 


+  a(^+^+1  ,+  ,  .  -  ^+J;  , ) 

i,j+l,k  i,j,k+l  i,j,ky 


,n+l 


,.n+l 


(13a) 


Vn+1. 


TIT  I  *  A  X  JT 

i.j-l.k-1  =  *1, j-l,k-l'  a(ti,j-l,k  +  Vj,k-l'  'J'i,j,k) 


+  ot('p”+1.  ,  .  +  ^  .  -  V**1.  .  ) 

i , J -1 . k  i,j,k-l  i , J , k 


,n+l 


,.n+l 


(13b) 


^  'p  «Hf 

where  ¥.  .  ,  =  [U,  .  ,  ,  U,  ,0]  with  U  and  W  determined  from 

1»J,K  ltJ,K  1 »  J . k 

previous  marching  locations  and  0  £  a  £  1.  The  consistent  form  of  the 
solution  algorithm  is  obtained  when  of=1.  After  substitution  the  terms 
in  equation  (11)  may  be  regrouped  giving  relations  for  GM .  ,  ,  ET  .  ,  ,  and 

J  ,  K  J  ,  K 

FT  .  ,  .  A  simple  two  pass  inversion  procedure  may  now  be  used  to 

j ,  K 

determine  the  solution  in  the  cross  plane. 

The  CCSIP  in  its  present  form,  equation  (12),  is  not  a  symmetric 
algorithm.  This  asymmetry  in  the  algorithm  can  generate  asymmetries  in 
the  cross  plane  solution  for  three-dimensional  geometries.  These 
asymmetries  in  the  cross  plane  solution  can  lead  to  a  procedure  that  is 
globally  divergent.  Symmetric  solutions  may  be  restored  by  iteration  on 
the  corner  points  (13a,b).  A  large  number  of  iterations  at  each  cross 
plane  is  required  to  restore  the  symmetry  of  the  solution,  however.  An 
alternative  means  of  obtaining  symmetric  solutions  is  to  make  the  CCSIP 
algorithm  symmetric  by  solving  the  mirror  image  of  equation  (12).  The 
mirror  image  solution  may  then  be  averaged  with  the  solution  from 
equation  (12).  This  averaging  produces  a  symmetric  solution  in  the 
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cross  plane  while  requiring  work  corresponding  to  only  one  local 
iteration. 

5.3  Initial  Conditions  and  Convergence  Acceleration 

We  have  found  the  present  relaxation  procedure  to  be  sensitive  to 
the  initial  guess  for  4>;  the  assumption  of  uniform  flow  is  generally 
inadequate  as  an  initial  condition.  Problems  associated  with  the 
initial  conditions  may  be  remedied  by  starting  the  viscous  solution  with 
the  corresponding  solution  to  the  inviscid  full  potential  equation. 
This  initial  guess  may  be  obtained  by  solving  the  inviscid  form  of  the 
composite  velocity  equations  with  slip  boundary  conditions.  An 
attractive  alternative  to  using  the  composite  system  for  the  potential 
flow  is  to  use  a  more  efficient  code  developed  strictly  for  the  full 
potential  equation.  The  composite  velocity  procedure  may  then  be  used 
to  efficiently  determine  the  viscous  solution  to  the  problem. 

The  convergence  of  the  global  relaxation  procedure  for  the  composite 
velocity  system  is  determined  primarily  by  the  convergence  rate  for  the 
inviscid  flow  variable,  <J>.  It  is  well  known  that  plane  relaxation 
procedures  slow  down  dramatically  as  the  number  of  planes  in  the 
marching  direction  is  increased.  In  the  present  calculations,  solutions 
are  always  obtained  on  the  coarsest  possible  grid;  subsequent  solutions 
on  finer  grids  are  then  initiated  with  the  coarser  grid  solutions.  This 
allows  for  more  rapid  upstream  propagation  of  downstream  information. 

The  previous  discussion  suggests  that  convergence  could  be  greatly 
accelerated  using  a  formal  multi -grid  procedure.  Multi -grid  procedures 
are  particularly  appealing  for  the  composite  velocity  system  since  the 
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inviscid  solution  is  based  on  the  full  potential  equation.  A  uni¬ 
directional  multi -grid  procedure  has  already  been  developed  by  Himansu 
and  Rubin[15]  for  the  RNS  equations  in  primitive  variable  form.  This 
technique  can  easily  be  adapted  for  the  composite  velocity  system  and  is 
recommended  for  further  investigation. 

The  convergence  of  the  potential  flow  may  also  be  accelerated  by 
using  overrelaxation  for  the  calculation  of  A4>  in  the  continuity  (i.e., 
full  potential)  equation  [16].  This  introduces  a  time  like  term  <t> 

s t 

that  provides  temporal  damping  to  the  system.  This  additional  damping 
aids  in  accelerating  the  convergence  of  the  inviscid  flow  while 
introducing  no  additional  error  in  the  converged  state. 

5.4  Transonic  Flow  and  Turbulence  Model 


Since  the  outer  inviscid  solution  of  the  composite  velocity  system 
is  based  on  the  full  potential  equation,  techniques  for  treating  the 
transonic  full  potential  equation  may  be  adopted  for  the  composite 
velocity  system.  The  flux  biasing  procedure  of  Enquist-Osher [7] ,  which 
specifies  a  modified  density  for  the  full  potential  equation,  is  used 
for  the  present  calculations.  This  is  a  simple  scheme  that  has  several 
advantages.  The  scheme  provides  monotone  profiles  through  captured 
shock  waves.  The  captured  shocks  are  very  sharp  with  only  a  two  point 
transition  through  the  shock  wave.  The  flux  biasing  scheme  also 
accurately  monitors  sonic  conditions  while  transitioning  through  sonic 
lines.  The  theoretical  reason  for  these  properties  lies  in  the  fact 
that  the  Enquist-Osher  scheme  satisfies  an  "entropy  inequality." 
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The  numerical  solution  of  the  RNS  equations  for  large  Reynolds 
number  flows  requires  the  introduction  of  a  model  for  turbulence. 
Algebraic  eddy  viscosity  models  provide  the  simplest  means  for  modelling 
turbulent  flows.  Implementation  of  these  models  requires  minimal 
computer  time  and  storage.  The  algebraic  eddy  viscosity  model  used  in 
the  present  work  was  developed  by  Baldwin  and  Lomax[17j.  This  two-layer 
model  is  patterned  after  the  Cebeci -Smith  turbulence  model,  but  removes 
the  need  for  determining  the  edge  of  the  boundary  layer.  For  the 
afterbody  configurations  calculated  herein,  the  Shang  and  Hankey 
relaxation  model  [18]  has  been  used  to  account  for  the  "memory"  of 
turbulence  transport.  This  is  required  for  flows  with  large  adverse 
pressure  gradients. 
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6.  Results  and  Discussion 


The  composite  velocity  procedure  will  be  used  to  obtain  results  for 
nonaxisymmetric  afterbody  configurations.  To  test  various  aspects  of 
the  code,  several  simpler  configurations  have  also  been  calculated, 
★ 

figure  1.  A  sampling  of  these  results  will  be  presented  here.  A  more 
extensive  discussion  of  these  results  is  given  in  references  20  and  21. 

6.1  Axisymmetric  Circular  Arc  Bump 

The  flow  over  an  axisymmetric  circular  arc  bump  is  used  to  test  the 
ability  of  the  present  technique  to  calculate  turbulent,  transonic 
flows.  The  particular  geometry  chosen  has  been  used  by  Johnson  [19]  to 
evaluate  the  capabilities  of  various  turbulence  models  to  predict 
inviscid-viscous  interactions  that  occur  at  transonic  conditions.  The 
test  model  consists  of  a  circular  cylinder  having  a  diameter  of  15.2  cm 
with  a  circular  arc  bump  having  a  maximum  thickness  of  1.9  cm  and  a 
chord  length  of  20.3  cm.  To  eliminate  problems  from  grid 
discontinuities  at  the  leading  and  trailing  edge  of  the  bump,  the  sharp 
corners  have  been  smoothed  using  a  third  order  polynomial. 

In  the  axial  direction,  a  uniform  spacing  A|=0. 01667  is  specified  on 
the  circular  arc  bump.  A  uniform  stretching  factor  of  o=1.2  is 
specified  ahead  and  aft  of  the  bump  with  the  grid  extending  4  chord 
lengths  in  each  direction  and  a  maximum  grid  spacing  of  A£=0.7.  A  total 


★ 

Figures  are  located  at  end  of  report. 
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of  105  grid  points  is  used  in  the  axial  direction  with  61  points  on 
the  bump  and  22  grid  points  ahead  and  aft  of  the  bump.  An  initial 
spacing  at  the  wall  of  At}=0. 000025  is  used.  A  total  of  61  grid  points 
is  used  in  the  normal  direction.  The  grid  spacing  at  the  wall  gives  a 

value  of  rj+ ~  1  for  the  first  point  away  from  the  wall.  The  outer 
boundary  is  located  6.2  radii  from  the  surface  of  the  cylinder. 

Results  for  a  supercritical  case,  M  =0.85  and  Re=2. 03x10*  based  on 

chord  length  are  given  in  figures  2  and  3.  Comparison  of  the  pressure 
coefficient  with  the  experimental  results  of  Johnson  shows  excellent 
agreement  of  the  shock  strength  and  location.  Downstream  of  the  shock, 
the  calculated  pressure  shows  a  larger  compression  than  the  experimental 
results  in  the  trailing  edge  region.  From  the  wall  shear  results  this 
is  a  region  of  reverse  flow.  The  computed  separation  region  lies 
between  x=0.86  and  x=1.04.  The  experimental  results  show  a  separated 
region  of  flow  between  x=0. 815-0. 85  and  x=l.l.  The  variance  in  the 
separation  point  in  the  experimental  results  arises  from  the  fact  that  a 
rapid  increase  in  the  size  of  the  separation  bubble  occurs  at  M^O.85. 

The  experimental  results  show  a  slightly  larger  separated  region  than 
the  computed  results.  The  discrepancy  in  the  size  of  the  separated  flow 
region  and  the  pressure  results  in  this  region  is  due  uo  the  inability 
of  the  turbulence  model  to  adequately  treat  separated  flow  regions.  The 
smoothing  of  the  trailing  edge  may  also  be  causing  the  computed  flow  to 
rr attach  earlier  than  the  experimental  results. 
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6.2  Elliptical  Cross  Section  with  Circular  Arc  Axial  Variation 

For  the  first  three-dimensional  flow  calculations,  a  cylinder  with  a 
2:1  elliptical  cross  section  is  used.  A  circular  arc  variation  in  the 
axial  direction  is  specified  with  the  major  axis  having  the  same  axial 
variation  as  the  previous  axisymmetric  geometry.  A  grid  size  of 
109x61x17  is  used  to  calculate  the  present  results.  A  uniform  spacing 
of  A£=0. 01667  is  specified  on  the  bump  with  a  grid  stretching  factor  of 

1.2  specified  ahead  and  aft  of  the  bump.  The  grid  is  stretched  to  a 
maximum  grid  spacing  of  Af=0.125  which  is  then  held  constant.  The  grid 
extends  2.27  chord  lengths  ahead  and  aft  of  the  body.  In  the  azimuthal 
direction,  a  uniform  spacing  A<f=0.098  radians  is  specified.  The  normal 
grid  distribution  is  the  same  as  that  used  for  the  axisymmetric 
calculations. 

The  flow  calculated  is  for  Re=2. 03x10*  and  M  =0.9.  A  turbulent 

oo 

boundary  layer  profile  is  specified  at  the  inflow.  The  axial  variation 
of  the  wall  shear  and  the  pressure  coefficient  for  several  azimuthal 
locations  is  given  in  figures  4  and  5.  The  wall  shear  shows  a 
substantial  variation  in  the  azimuthal  direction  while  the  pressure 
coefficient  shows  less  variation  except  near  the  leading  and  trailing 
edge.  There  is  a  sharp  decrease  in  wall  shear  around  x=0.3  and  then  a 
plateau  region  before  the  decrease  in  wall  shear  owing  to  the  shock. 
This  behavior  is  due  to  the  curvature  effects  of  the  geometry  on  the 
boundary  layer  near  this  symmetry  plane.  The  shock  strength  and 
location  show  very  little  variation  in  the  azimuthal  direction. 
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The  effects  of  grid  resolution  in  both  the  axial  and  azimuthal 
direction  were  also  investigated  for  this  geometric  configuration.  The 
predominant  effect  of  grid  refinement  is  on  the  wall  shear  distribution 
in  the  region  of  the  shock.  Grid  refinement  has  a  minimal  effect  on  the 
pressure  coefficient.  The  effect  of  second  order  accuracy  for  the 

and  derivatives  has  also  been  investigated  with  the  main  effect  again 


being  on  the  wall  shear  distribution  in  the  region  of  the  shock.  A  more 
detailed  discussion  of  these  results  is  again  found  in  references  [20- 
21]- 


6.3  Nonaxlsymmetrlc  Afterbody  Configuration 

The  nonaxi symmetric  afterbody  configuration  for  which  calculations 
are  to  be  made  is  shown  in  figure  6.  The  geometry  consists  of  a  main 

body  section  with  an  elliptical  cross  section  of  aspect  ratio  1.097. 

The  upper  symmetry  plane  transitions  to  a  final  boattail  angle  of  17.6°, 
while  the  side  symmetry  plane  transitions  to  a  final  boattail  angle  of 
0 

6.7  .  A  solid  sting  is  used  to  simulate  the  nozzle  exhaust  plume.  The 

final  aspect  ratio  of  the  elliptical  cross  section  is  1.94.  A  third 

order  polynomial  is  used  to  smooth  the  discontinuity  at  the  juncture  of 
the  boattail  with  the  sting. 

The  flow  fields  calculated  are  doubly  symmetric  and  only  one  quarter 
of  the  flow  field  is  calculated.  A  grid  of  146x65x17  grid  points  is 
used.  Twenty-eight  planes  are  used  to  resolve  the  boattail  with  the 
spacing  on  the  boattail  being  0.4%  of  the  total  body  length.  Sixty-five 
grid  points  are  used  in  the  body  normal  direction  with  a  minimum  spacing 
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at  the  wall  of  0.001%  of  the  body  length  and  the  outer  boundary  located 
approximately  one  body  length  away. 

The  first  case  calculated  for  this  geometry  is  for  M  =  0.8  and 

6 

Re=1.0  x  10  .  The  flow  remains  completely  subsonic  for  this  freestream 
Mach  number.  The  pressure  coefficient  and  wall  shear  distribution  for 
various  rj=c onst  lines  are  given  in  figures  7a, b.  The  beginning  of  the 
boattail  is  located  at  x=0.88.  Large  azimuthal  variations  in  the 
pressure  and  wall  shear  on  the  boattail  are  seen.  At  the  point  of 
minimum  pressure,  there  is  a  20%  variation  in  pressure  from  the  top  to 
side  symmetry  plane.  The  point  of  minimum  pressure  on  the  side  symmetry 
plane  occurs  slightly  forward  of  the  minimum  pressure  location  on  the 
top  symmetry  plane. 

The  methods  for  accelerating  the  calculations  described  earlier  were 
used  for  these  results.  A  total  of  1  hour,  25  minutes,  and  22  seconds 
of  CPU  time  on  the  Cray  XMP  were  required  for  this  calculation.  Of  the 
total  time,  52  minutes  and  11  seconds  or  61%  of  the  time  is  spent 
calculating  the  potential  flow  used  as  an  initial  condition  for  the 
viscous  calculation.  This  time  can  be  greatly  reduced  by  using  a  code 
designed  specifically  for  calculating  potential  flows.  Thus  with  a 
potential  flow  solution  to  the  problem,  a  viscous  flow  solution  may  be 
obtained  in  only  33  minutes  and  11  seconds.  The  amount  of  memory 
required  for  this  calculation  is  907,776  words. 

The  pressure  coefficient  and  wall  shear  distribution  for  a  transonic 

case,  Mm=  0.9,  are  given  in  figures  8a, b.  A  large  azimuthal  variation 

in  both  the  pressure  and  the  wall  shear  distribution  can  again  be  seen. 
The  shock  wave  on  the  top  symmetry  plane  is  much  stronger  and  is  located 
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further  aft  on  the  boattail.  A  24%  variation  in  the  minimum  pressure  is 
seen.  A  small  region  of  reverse  flow  induced  by  the  shock  wave/boundary 
layer  interaction  can  be  seen  in  the  wall  shear  distribution  just 
downstream  of  the  shock  wave. 

Plots  of  pressure  and  Mach  contours  show  a  clearer  picture  of  the 
flowfield  characteristics.  Figure  9  gives  surface  pressure  contours  in 
the  boattail  region.  The  expansion  region  terminating  in  a  shock  wave 
can  clearly  be  seen.  This  region  has  a  greater  axial  extent  along  the 
top  symmetry  plane,  and  the  shock  wave  again  appears  much  stronger. 
Pressure  contours  in  two  cross  planes  located  at  x=0.909  and  x=0.934  are 
given  in  figures  10  and  11  respectively.  These  again  show  the  three- 
dimensional  effects  in  the  flowfield.  In  particular,  there  is  a  very 
large  azimuthal  variation  in  the  pressure  in  the  x=0.934  cross  plane, 
which  is  located  in  the  shock  region. 

Mach  contours  in  the  top  and  side  symmetry  planes  are  given  in 
figures  12  and  13  respectively  which  clearly  show  the  sharp  shock  waves 
generated  by  the  Enquist-Osher  flux  biasing  procedure.  The  shock  is 
again  much  stronger  on  the  top  symmetry  plane.  The  effect  of  the  shock 
on  the  boundary  layer  can  also  be  seen  with  the  boundary  layer  becoming 
thinner  in  the  region  ahead  of  the  shock  and  thickening  as  the  boundary 
layer  interacts  with  the  shock  wave. 

Figures  14  and  15  show  the  effect  of  viscosity  on  the  pressure 
coefficient  along  the  top  and  side  symmetry  planes.  The  viscous  effects 
tend  to  weaken  and  diffuse  the  shock  with  the  viscous  shock  being 
located  slightly  upstream  of  the  potential  shock.  The  viscous  flow  also 
undergoes  a  much  smaller  compression  in  the  boattail  juncture  region. 
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6 ■  Future  Work 


We  have  noted  several  possible  needed  improvements  to  the  current 
composite  velocity  scheme  during  the  course  of  these  calculations.  An 
improved  downstream  boundary  condition  that  more  accurately  represents 
the  three-dimensional  effects  at  this  boundary  is  required.  We  have 
found  the  solution  procedure  to  be  sensitive  to  the  downstream  boundary 
condition  in  a  very  local  region  near  this  boundary.  The  current 
implementation  of  the  v=0  boundary  condition  on  the  body  also  needs  to 
be  improved  so  that  v=0  is  maintained  for  each  global  iteration.  This 
should  improve  the  convergence  performance  of  the  algorithm.  Finally, 
work  is  currently  under  way  to  replace  the  CCSIP  algorithm  with  a  direct 
sparse  matrix  solver.  This  should  further  improve  the  robustness  of  the 
present  numerical  scheme. 
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Nomenclature 


C 

P 

c 

P 

Fkleb 

1  j 

Sij.  S 
s 

G 

H 

M 

P 

q 

Re 

S 

SM.SP 

T 

u 

U 


u 

i 

V 

w 

W 


skin  friction  coefficient 
pressure  coefficient 
specific  heat  at  constant  pressure 
Klebanoff  intermittency  factor 


metric  tensors 


determinant  of  the  metric  tensor 


Bernoulli  pressure  term 

total  enthalpy 

Mach  number 

pressure 

speed 

Reynolds  number 
entropy 

switches  for  forward  and  backward  differencing 
temperature 

velocity  component  in  axial  direction 

viscous  composite  velocity  component  in  axial 

direction 

inviscid  velocity  in  axial  direction 


velocity  component  in  body  normal  direction 
velocity  component  in  azimuthal  direction 
viscous  composite  velocity  component  in  azimuthal 
direction 
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V 

e 


x.y.z 


y 


e 


t'V't 


inviscid  velocity  in  azimuthal  direction 
Cartesian  coordinates 

factors  for  stretched  grid  in  normal  direction 

ratio  of  specific  heats 

Christoff el  symbol  of  the  second  kind 

switch  for  first  or  second  order  accuracy 
coordinates  in  computational  plane 


+ 

V 


law  of  the  wall  coordinate 


0 


H,\ 


P 

a 

* 

A<t> 

* 


azimuthal  location  -  Jy*  +  z* 
coefficients  of  viscosity 
eddy  viscosity 


density 

grid  stretching  factor 
potential  like  function 


<j>n+1  .  <j>n 


[U.  W,  A*]' 


vorticity 
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Subscripts 


t,  x, y, z 

i .  j .  k 


.  j 


OO,  f  S 

m 

w 


derivatives 

indexes  for  coordinate  directions 

or  counters  for  tensor  notation 

covariant  derivative 

freestream 

maximum 

wall  value 


Superscripts 


i.  j,k 

n+1 

n 


counters  for  tensor  notation 
new  iteration  level 
old  iteration  level 
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m mm 

mn 


HkI 

BMmf 

Inal 
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